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An investigation is conducted into the determination of the credibility of 
interacting boundary layers in predicting compressible subsonic flows over 
smooth surface imperfections. The case of smooth backward-facing steps is 
considered. The predicted mean flows are compared with those obtained 
using a Navier-Stokes solver. Moreover, the linear 2-D compressible stability 
characteristics of both mean flows are compared. The results show that the 
interacting boundary-layer formulation produces accurate mean flows that 
yield accurate linear stability characteristics, such as growth rates and 
amplification factors. 
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I. Introduction 


The boundary layers over natural laminar flow components in the 
presence of surface imperfections (e.g., waviness and steps) must accurately 
be computed so that the effect of these imperfections on the stability and 
transition can be evaluated. Moreover, the magnitudes of the imperfections 
under consideration are such that strong viscid-inviscid interaction and small 
separation bubbles are unavoidable. Definitely, solutions to the full 
Navier-Stokes equations can accurately predict such flowfields provided that 
the grid is fine enough so that important flow structures are not smeared by 
the truncation errors or artificial dissipation. However, the number of flow 
cases that needs to be investigated is very large, and hence solving the full 
Navier-Stokes equations is a very expensive task. A more economical 
alternative is to solve the interacting boundary-layers equations. 

The purpose of the present paper is to investigate the accuracy of the 
compressible interacting boundary-layer computations and to establish their 
credibility in accurately predicting flows over surface imperfections. 
Specifically, we consider the flow past a flat plate having a smooth 
backward-facing step. Our approach is to compare the results of interacting 
boundary-layer computations with solutions to the Navier-Stokes equations. 
Comparisons are made for the mean-flow profiles as well as the stability 
characteristics, such as, the growth rates and amplification factors of linear 
stability waves. 


II. Mean Flow 

A. Navier-Stokes Solver 

The thin-layer compressible Navier-Stokes equations were solved using 
the well known computer code "ARC2D" which was developed at NASA Ames; 
the version that we acquired from NASA Ames is designated 1.5 GAMMA, 
7/2/85. The code incorporates different methods of solution, all of which are 
implicit in time, and uses second-order central differences in space. We 
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selected the method of solution in which the diagonal form of the equations is 
used'. Mixed second- and fourth-order dissipation terms were added 
explicitly and implicitly, and the resulting pentadiagonal system of equations 
were solved directly. 

Velocities are normalized using the freestream velocity U‘ x , lengths are 
normalized using the distance L’ from the leading edge of the plate to the 
center of the step, and the temperature, viscosity, and thermal conductivity 
coefficients are normalized using their freestream values 7^,^, and k^, 
respectively. Sheared Cartesian grids are used for all the cases analyzed in 
this paper. An example is shown in Fig. 1 for a smooth backward-facing step. 
The equation of the step is 


y = yM 1+erfC). (1) 

where 

£ = fte -3/8 /l 5/4 (x — 1), (2) 

Re is the Reynolds number based on the distance from the leading edge to the 
step center (x = x*jL* = 1), and X = 0.332057. The step specified by Eq. (1) was 
originally employed by Smith and Merkin 2 who analyzed the incompressible 
mean flow using triple-deck theory. The numerical values used in Fig. 1 are 
Re = 10® and h — — 0.003. We note that the y-coordinates in Fig. 1 are 
magnified by a factor of 20 relative to the x-coordinates. 

The inflow boundary of the computational domain is located at x^ — 0.06 
(i.e., the plate leading edge x = 0 is included in the domain) and the outflow 
boundary is located at x = 2.0. The top boundary is placed at yc;0.4. More 
details about the grid are given in Section IV. 

B. Interacting Boundary Layers 

We developed a code for solving the compressible interacting 
boundary-layer equations. The numerical method is an extension of 
Veldman's method 3 to compressible flows. The salient feature of the method 
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is the simultaneous solution of the boundary-layer equations and the inviscid 
flow, which is given by the small disturbance theory of compressible potential 
flow. A similar treatment was presented by Davis 4 for subsonic and supersonic 
flows over parabolic humps. For more details about our method we refer the 
reader to Ref. 5. 

The governing equations are the compressible boundary-layer equations. 
In terms of a combination of the Levy-Lees variables and the Prandtl 
transposition theorem, the problem is given by 

2(FF( + VI :^J-(e^ + ^-Q) = 0 (3) 

2£F { + V„ + F = 0 (4) 

2 CPQi + w, - -L ( ±r -S. ) - ft, - 1)Mi £. = 0 (5) 

subject to the boundary conditions 


F = 0, V = 0, Q„ = 0 at rj = 0 

( 6 ) 

F — ► 1 and Q — ► 1 as tj -* oo 

( 7 ) 

F(Z 0 , rj) and Q = Q(£ 0 , rj) at { = £ 0 

(8) 


where 


m = f 


P e P e U e dx and '/M = 




\ Z pdx 


r , ) = 


F = = = (x)] 


Pe^ePe 


( 9 ) 


( 10 ) 


( 11 ) 
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•*-f :^t (12 > 


and C* and c; are the coefficients of the specific heat of the gas at constant 
pressure and volume, respectively. Here p e , p 9 , and U 9 are the edge density, 
viscosity coefficient, and velocity, respectively. The displacement thickness is 
given by 


5 = 


1 


r- 

Pe^e *o 


(Q - F)drt 


(13) 


The interaction law relates the edge velocity U e to the displacement surface. 
Using thin airfoil theory, we obtain 


U e = U e + 


1 f°° 

pn J L£ x - 


<5 d(lnp e ) 


t 


dt 


+ 


1 

pn J Ll 


d(U e S) 


x - t 


dt 


dt 


(14) 


where 

/* = V 1 - < i5 > 

and U 9 is the inviscid surface velocity in the absence of the boundary layer, 
which is also determined using thin airfoil theory. After some manipulations, 
the interaction law can be rewritten as 


U e 


= 1 + 



1 1 r°° 

x-t dt + fin J L£ 


U e § 
x - 1 


<XJnPe) 

dt 


dt 


(16) 


where 


X = f+U e 8 (17) 

and the principal values of the Cauchy integrals in Eqs. (14) and (16) are 
assumed. Equations (3)-(8), (13), and (16) are solved simultaneously following 
the procedure of Veldman 3 . Veldman integrated Eq. (16) by parts to obtain a 
second derivative for x and expressed U 9 as a linear combination of the values 
of 5 at the nodes. However, we follow Davis and Werle 6 and Nayfeh et al 5 and 
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perform the integration by parts to eliminate the derivative of x and assume x 
to vary linearly over the differencing intervals, resulting in a second-order 
accurate scheme. The finite-difference scheme used employed a three-point 
upwind differencing formula for the £ derivatives and a central differencing for 
the rj derivatives. Upwinding the £ derivatives stabilizes the numerical method 
in regions of reversed flow (i.e., F<0). The system of equations is solved by 
a scheme similar to that employed by Nayfeh et al 5 . The distribution of the 
mesh points on the steps was chosen in accordance with the triple-deck 
scaling. The lower deck was resolved using a mesh with a variable step size 
in the 77 direction so that a larger number of points could be employed near the 
wall than in the rest of the boundary layer. 


III. Stability Formulation 

We consider the linear quasiparallel two-dimensional stability of 
two-dimensional compressible flows over smooth backward-facing steps. The 
quasiparallel assumption can be justified a posteriori 5 . The calculated 
wavelengths are the order of the boundary-layer thickness. We superpose a 
small time dependent disturbance on each mean-flow, thermodynamic, and 
transport quantity. Thus, we let 

q(x,y,t) = q m (y) + g(x,y,f) (18) 

where q m (y) is a two-dimensional basic-state quantity and q(x,y,t ) is a 
two-dimensional unsteady disturbance quantity. Here, q stands for the 
velocity components (u and v), temperature T, pressure p, density p, and 
viscosity coefficient p. Substituting Eq. (18) into the Navier-Stokes equations, 
subtracting the basic-state quantities, and linearizing the resulting equations 
in the q's, we obtain 


dp 


+ P 


du 

m dx 


+ Um dx + dy ~ ° 


dt 


(19) 



m = XJn m and r = 2 + m (25) 

The problem is completed by the specification of the boundary conditions; they 


u = v = T = 0 at y = 0 (26) 

u, v, T -> 0 as y -> oo (27) 

Next, we assume normal-mode solutions of the form 
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where q stands for the eigenfunctions of u, v, p, and T. Using Eq. (28) in Eqs. 
(19)-(23), (26), and (27) and dropping the tilde from q, we obtain the eigenvalue 
problem 


DT m 

Dv = — iau H — - — v 4- 


' m 


i£p jQj 

Pm T m 


(29) 


d 2 u = ( + m A u _ pu 


+ 


Pm 

p m RDu m p m 'DT t 


Pm 


P 


— IOL 


m 


m 


Pm 


v - /(I + m)ccDv 


^ p 




D (Pm') + 


D 2 u 


Pm 


m f 


T-^Du m DT 


(30) 


Xo Op= — ia(r-^~~+ 2y " PT ” 1 \u-iaDu + 


' m 


f^m 


iRQ. 2 
— a 


+ r OJ^ + r^prj \ + ._r_ 


' m 


Pm“^i 


m 


' m 


Pm^m 

DT, 


n 


m 
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' m 


Pm 

Pm 


DT, 


m 


1(31) 


- « Du m ~\p 4 - 


iaDuJ 4 s - + + 


' m 


i' ft y or 

Pm T m m 


t-J^Ldt 


' m 


-2(y- 

- 1 )M 2 PrDu m Du 4- r — - 

H-m 

L 

2 i(v - 

1)a/v£Pr0u m 

V 

+ i(y - 


PrR -& 

p+[- 

• /PrP n -4- 

r-m 

4- a 2 - 

(07 ) 2 
Pm { m) 

(32) 

Pm' 

f^m 

£> 2 ^ - (y - 

1 )MlPr 

Pm (n , 2 ‘ 
Pm 

7-2 

Pm DT m 
Pm 




u = 

v=T = 

0 at y = 0 



(33) 



a, u 

i 

o 

as y — ♦ co 



(34) 


where 



(35) 


D. = a> — au m 


andp ™= 


PjvTm 

vmL 


t*m = 


c/r, 


m 


= 




m 


dT, 


m 


(36) 


(37) 


Equations (29)-(32) can be written as a system of six first-order equations of the 
form 


DC = A( (38) 

where 

C T = {u Du v p T DT} t (39) 

The boundary conditions (33) and (34) can be written as 

Cl = C 3 = C 5 = 0 at y = 0 (40) 

C n -* 0 as y -+ oo (41) 

The system of equations (38) subject to the boundary conditions (40) and 
(41) constitutes an eigenvalue problem. For a given Reynolds number, Mach 
number, and mean-flow profiles, we determine the dispersion relation 

g(co, a) = 0 (42) 

employing a numerical procedure. For the spatial stability problem 
considered here, we specify co and an initial guess for the eigenvalue a. We 
numerically integrate the system (38) from y = y e to y = 0, where y„ indicates 
a value of y outside the boundary layer. In performing these calculations we 
used the computer code 'SUPORT' developed by Scott and Watts 7 . This code, 
which is based on the method of Gudonov, solves stiff two-point 
boundary-value problems. It uses the Runge-Kutta-Fehlburg scheme to 
integrate the equations and the Gram-Schmidt orthonormalization scheme to 
keep the solution vectors linearly independent. Normally, the first guess of the 
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eigenvalue is incorrect and therefore one of the boundary conditions at y = 0 
is not satisfied. Then, a Newton-Raphson scheme is used to iterate on this 
unsatisfied boundary condition and obtain the eigenvalue to the desired 
accuracy. 

Having found the growth rates — a„ we compute the N-factor as 


N = - 


'R 

2ctjdR 

J *o 


where R 0 corresponds to the Reynolds number at the first neutral point. 
During the course of integration, if we find that N becomes negative we set it 
equal to zero. 


IV. Comparison Between Interacting Boundary Layers and 
Navier-Stokes Results 

In this section solutions for flows over smooth backward-facing steps 
obtained using the Navier-Stokes solver and 1BL computations are compared. 
The Mach numbers considered for the comparison are M = 0.5 and M = 0.8. 
Comparisons are made for the mean flows as well as their stability 
characteristics, such as the growth rates and the N-factors. To predict the 
flowfield and its stability accurately, we must choose a proper grid, which is 
fine enough to capture the important flow structures. The far field boundary 
also plays an important role in the determination of the solution. The 
boundary must be far enough to prevent contamination of the solution by 
reflection there, especially in the high subsonic Mach number case. 

The grid used in the IBL calculations has a uniform step size of 
Ax = 0.005 in the streamwise direction and a geometrically stretched grid in 
the rj direction . At the wall A >7 = 0.05 and the stretching factor is 1.05. The 
edge of the boundary layer is taken at rj e = 8.0 . 
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Four grids were used with the Navier-Stokes solver. In grid 1 (136x70), a 
uniform Ax = 0.005 was used in the range 0.9<x<1.1, where strong 
streamwise gradients are expected. For x<0.9 and x>1.1,Ax is stretched 
geometrically at the rate of 1.05 provided that Ax did not exceed 0.03. If Ax 
exceeded 0.03, a uniform spacing Ax = 0.03 was used. To cover the domain 
—0.06 < x < 2.0, 136 points were used. The step size in the y direction was 
geometrically stretched between the wall and the top boundary y»0.4 with 
Ay, = 1.5x10-*, and the stretching factor was adjusted so that the number of 
points was 70. In grid 2 (136,99), the streamwise distribution of grid points was 
the same as in grid 1, while the y direction had 99 points with Ay, = 0.3x10-". 
In grid 3 (166x120), the streamwise distribution in the interval 0.9 <, x < 2.0 was 
the same as that in grid 1, but more points were added in the range 
-0.06 < x < 0.9 so that the total number of x-points was 166. The purpose of 
this addition is to obtain a finer grid near the leading edge of the plate. The 
y distribution was also refined so that 120 points were used between the wall 
and the top boundary y%0.4 with Ay, = 0.3x10 *. Grid 4 was the finest grid 
used, where 176 points were used in the x-direction and 153 in the y-direction. 
In this grid, the streamwise distribution of grid points was the same as that in 
grid 1 except that we started with Ax = 0.003 instead of 0.005. The y grid had 
153 points between the wall and the top boundary, which is placed now at 
y~0.8 instead of the 0.4, which was used for the first three grids, and 
Ay, = 0.3x10-" . 

First we compare the mean-flow characteristics predicted by the IBL code 
with those obtained by the Navier-Stokes code. The friction-coefficient 
distributions on a backward-facing step with h = - 0.003, Re = 10®, and 
M x = 0.5 are shown in Fig. 2. The Navier-Stokes results were obtained using 
grids 1, 2, and 3. The corresponding pressure-coefficient distributions are 
shown in Fig. 3. We note that there is a great discrepancy between the IBL 
computations and the Navier-Stokes results when grid 1 (corse grid) is used. 
The discrepancy is obvious downstream of the separation point. That 
prompted grid refinements in the Navier-Stokes calculations. The results of 
the IBL computations are in good agreement with the results obtained using 
grids 2 and 3 almost everywhere except in the reattachment region. However, 




the excellent agreement between the results obtained using grids 2 and 3, with 
grid 3 being the finer of the two, suggested that a further grid refinement may 
not result in an appreciable change in the reattachment region. And hence 
we proceeded with the stability analysis. 

We fixed the frequency <y at a value of 80, and solved the eigenvalue 
problem for the wavenumber a. The mean flow was obtained by using the IBL 
code or by using the ARC2D code using grids t, 2, and 3. The growth rates 
( - a,) are shown in Fig. 4 as a function of R = JWe ~ for the four sources of the 
mean flow. The corresponding N-factors are shown in Fig. 5. We note that the 
IBL results are far off from the results of grid 1 both for the growth rate and the 
N-factor. A dramatic improvement was obtained when the finer grids 2 and 3 
were used. Figure 5 shows that the IBL results are in better agreement with 
the results of grid 3 than with those of grid 2, especially upstream of the step. 
However, the IBL predicts N-factors in the reattachment region that are lower 
than those predicted by using the two grids. 

Next we decided to use a still finer grid in order to improve the agreement 
in the reattachment region. In grid 4 (176x153), we refined the grid notably in 
the x-direction, while the distribution in the y-direction was almost the same 
as in grid 3. In Fig. 6, we depict the friction coefficient distributions obtained 
by using the IBL code and the ARC2D code with grids 1, 2, 3, and 4. There is 
an obvious improvement in the agreement of the IBL and Navier-Stokes results 
in the reattachment region when grid 4 is used. In Fig. 7 we show the growth 
rates ( — a y ) obtained by using the IBL mean flow and the ARC2D code with 
grid 4. The N-factors are depicted in Fig. 8, they show considerable 
improvement in the agreement between the IBL computations and the 
Navier-Stokes computations using grid 4. 

A grid refinement study was also conducted at a Mach number of 0.8. The 
mean flow and stability results show trends similar to those for the case of 
Mach number of 0.5. The growth rates predicted by using the IBL and N-S 
mean flows are depicted in Fig. 9 for a backward-facing step of h = - 0.003, 
Re = 10 s , = 0.8, and F = 50x10 6 . The N-S mean flow is obtained by using 

grid 3. The corresponding N-factors are shown in Fig. 10. The agreement 
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between the two methods of obtaining the mean flow is good. For more 
details, we refer the reader to Ref. 8. 

We conclude that for the purpose of the stability analysis of boundary 
layers over smooth surface imperfections, which might induce small 
separation bubbles, the interacting boundary-layer formulation is a viable 
alternative to the Navier-Stokes equations. 

V, Effects of Mach Number 

Using interacting boundary-layer mean flow computations, we present in 
Figs. 11 and 12 the pressure coefficient and skin-friction coefficient, 
respectively, for a backward-facing step of h = - 0.003, Re = 10 6 , and 

= 0.0, 0.5, and 0.8. It is evident in Fig. 12 that increasing the Mach number 
increases the streamwise extent of the separation region. For each Mach 
number, we determined the most amplified frequency, that is, the frequency 
which results in the highest N-factor. For M x = 0.0, 0.5, and 0.8, the 
corresponding frequencies are (60, 55, and 50) x 10 6 , respectively. The 
growth rates and N-factors are shown in Figs. 13 and 14, respectively. As 
expected, the maximum growth rate is reduced by compressibility, but 
because of the increase in the separation region with increasing Mach 
number, the growth rate in the reattachment region is higher for high Mach 
numbers. Although the maximum N-factor is reduced by increasing the Mach 
number, the N-factor at the end of separation is practically equal for the three 
Mach numbers. The increase in the separation region offsets the stabilizing 
effects of compressibility. 
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Figure Captions 


Figure 1 . A typical computational grid. 

Figure 2. Influence of grids on the skin-friction coefficient for M ^ = 0.5 and Re 

= 1.0x10 s : 1 IBL, - - - 2 NS (GRID1), 3 NS (GRID2), and ... 

4 NS (GRID3). 


Figure 3. Influence of grids on the pressure coefficient for M x = 0.5 and 

Re = 1.0x10 s : 1 IBL, - - - 2 NS (GRID1), 3 NS (GRID2), and 

... 4 NS (GRID3). 

Figure 4. Influence of grids on the growth rates for = 0.5 at Re = 1.0x10 s 

and F = 80x10“® : 1 IBL, - - - 2 NS (GRID1), 3 NS (GRID2), 

and ... 4 NS (GRID3). 

Figure 5. Influence of grids on the N-factor for M x = 0.5 at Re = 1.0x10 s and F 

= 80x10-® : 1 IBL, 2 NS (GRID1), 3 NS (GRID2), and ... 

4 NS (GRID3). 

Figure 6. Comparison of the skin-friction coefficient calculated by using 
interacting boundary layers with those calculated by using the 
Navier-Stokes solver and different grid refinements for IW M = 0.5 and 

Re = 1.0x10 s : 1 IBL, - - - 2 NS (GRID1), 3 NS (GRID2), ... 4 

NS (GRID3), and 5 NS (GRID4). 

Figure 7. Comparison of the growth rates based on the profile calculated by 
using the interacting boundary layer code with those based on the 


is 


profile calculated by using the Navier-Stokes solver with the finest 

grid for = 0.5 , Re = 1.0x10 s , and F = 80x10 6 : 1 IBL and 

2 NS (GRID4). 

Figure 8. Comparison of the N-factors based on the profile calculated by using 
the interacting boundary layer code with those based on the profiles 
calculated by using the Navier-Stokes solver with different grid 

refinements for = 0.5 , Re = 1.0x10 s , and F = 80x10”® : 1 IBL, 

- - - 2 NS (GRID1), 3 NS (GR1D2), ... 4 NS (GRID3), and 

5 NS (GRID4). 

Figure 9. Comparison of the growth rates based on the profile calculated by 
using the interacting boundary layer code with those based on the 
profile calculated by using the Navier-Stokes solver with the finest 

grid for M x = 0.8 , Re = 1.0x10 s , and F = 50x10-® : IBL, and 

NS (GRID4). 

Figure 10. Comparison of the N-factors based on the profile calculated by 
using the interacting boundary layer code with those based on the 
profile calculated by using the Navier-Stokes solver with the finest 

grid for = 0.8 , Re = 1.0x10®, and F = 50x10” s : IBL and 

NS (GRID4). 

Figure 11. Influence of Mach number on the pressure coefficient at Re = 
1.0x10 s : /W = 0.0, - - - M = 0.5 , and M = 0.8. 

Figure 12. Influence of Mach number on the skin-friction coefficient at Re = 
1.0x10 s : M 00 = 0.0, - - - M„ = 0.5 , and = 0.8 . 
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Figure 13. Influence of Mach number and most dangerous frequency on the 

growth rates for Re = I.OxlO 6 : M x = 0.0 and F = 60x10- 6 ; 

IW M = 0.5 and F = 55x1 0 -6 ; and IW M = 0.8 and F 

= 50x1 0- 8 . 

Figure 14. Influence of Mach number and most dangerous frequency on the 

N-factors for Re = 1.0x10 s : M x = 0.0 and F = 60x1 0” 8 ; 

= 0.5 and F = 55x1 0 -8 ; and M x = 0.8 and F = 50x1 0 -6 . 
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